The Richardson's Law in Large-Eddy 
Simulations of Boundary Layer flows 

G. Gioia 1 , G. Lacorata 1 , E.P. Marques Filho 2 , 
A. Mazzino 1 ' 3 and U.Rizza 1 

1 ISAC-CNR, Sezione di Lecce, 1-73100, Lecce, Italy 
2 Institute of Astronomy, Geophysics and Atmospheric Sciences, 
University of Sao Paulo, 05508-900, Sao Paulo, Brasil 
3 Dipartimento di Fisica, Universita di Genova, 1-16146, Genova, Italy 

February 2, 2008 

Abstract 

Relative dispersion in a neutrally stratified planetary boundary 
layer (PBL) is investigated by means of Large-Eddy Simulations (LES). 
Despite the small extension of the inertial range of scales in the simu- 
lated PBL, our Lagrangian statistics turns out to be compatible with 
the Richardson t 3 law for the average of square particle separation. 
This emerges from the application of nonstandard methods of analy- 
sis through which a precise measure of the Richardson constant was 
also possible. Its values is estimated as Ci ~ 0.5 in close agreement 
with recent experiments and three-dimensional direct numerical sim- 
ulations. 

1 Introduction 

One of the most striking features of a turbulent planetary boundary layer 
(PBL) is the presence of a wide range of active length scales. They range 
from the smallest dynamically active scales of the order of millimeters (the 
so-called Kolmogorov scale), below which diffusive effects are dominant, to 
the largest scales of the order of ten kilometers. Such a large range of excited 
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scales are essentially a continuum and the distribution of energy scale-by- 
scale is controlled by the famous Kolmogorov's 1941 prediction (see Frisch, 
1995 for a modern presentation). 

One of the most powerful concepts which highlighted the dynamical role of 
the active scales in the atmosphere was due to Richardson (1926). He intro- 
duced in his pioneering work the concept of turbulent relative dispersion (see 
Sawford, 2001 for a recent review) with the aim of investigating the large 
variations of atmospheric turbulent diffusion when observed at different spa- 
tial scales. 

In his work, Richardson proposed a diffusion equation for the probability den- 
sity function, p(r, t), of pair separation. Assuming isotropy such an equation 
can be cast into the form 



dp(r,t) 1 d 
dt r 2 dr 



r D(r) — — — 



(1) 



where the scale-dependent eddy-diffusivity D{r) accounts for the enormous 
increase in observed values of the turbulent diffusivity in the atmosphere. 
The famous scaling law D[r) oc r 4//3 was obtained by Richardson (1926) 
from experimental data. From the expression of D(r) as function of r and 
exploiting Eq. (^J) the well known non-Gaussian distribution 

p(r, t) oc r 9/2 exp (-Cr 2/3 /t) (2) 

is easily obtained. 

This equation implies that the mean square particle separation grows as 

R 2 (t) = (r 2 (t)) = C 2 et 3 (3) 

which is the celebrated Richardson's "t 3 " law for the pair dispersion. Here 
C 2 is the so-called Richardson constant and e is the mean energy dissipation. 
Despite the fact that the Richardson's law has been proposed since a long 
time, there is still a large uncertainty on the value of C 2 . Some authors have 
found C 2 ranging from ~ 10~ 2 to ~ 10 _1 in kinematic simulations (see, for 
example, Elliot and Majda, 1996; Fung and Vassilicos, 1998), although for 
kinematic models an energy flux e can hardly be defined. On the other hand, 
a value C 2 ~ 0(1) (and even larger) follows from closure predictions (Monin 
and Yaglom, 1975). More recently, both an experimental investigation (Ott 
and Mann, 2000) and accurate three-dimensional direct numerical simula- 
tions (DNS) (Boffetta and Sokolov, 2002) give a strong support for the value 
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C 2 ~ 0.5. 

The main limitation of the state-of-the-art three-dimensional DNS is that 
the achieved Reynolds numbers are still far from those characterizing the 
so-called fully developed turbulence regime, that is the realm of the Richard- 
son's (1926) theory. Moreover, initial and boundary conditions assumed in 
the most advanced DNS are, however, quite idealized and do not match those 
characterizing a turbulent PBL, the main concern of the present paper. 
For all these reasons we have decided to focus our attention on Large-Eddy 
Simulations (LES) of a neutrally stratified PBL and address the issue related 
to the determination of the Richardson constant C 2 . The main advantage of 
this strategy is that it permits to achieve very high Reynolds numbers and, 
at the same time, it properly reproduces the dynamical features observed in 
the PBL. 

It is worth anticipating that the naive approach which should lead to the 
determination of C 2 by looking at the behavior of R 2 (t) versus the time t is 
extremely sensitive to the initial pair separations and thus gives estimations 
of the Richardson's constant which appear quite questionable (see Fig. 3). 
This is simply due to the fact that, in realistic situations like the one we 
consider, the inertial range of scales is quite narrow and, consequently, there 
is no room for a genuine t 3 regime to appear (see Boffetta et al., 2000 for 
general considerations on this important point). 

This fact motivated us to apply a recently established 'nonstandard' analysis 
technique (the so-called FSLE approach, Boffetta et al., 2000) to isolate a 
clear Richardson regime and thus to provide a reliable and systematic (that 
is independent from initial pair separations) measure for C 2 . This is the main 
aim of our paper. 

2 The LES strategy 

In a LES strategy the large scale motion (that is motion associated to the 
largest turbulent eddies) is explicitly solved while the smallest scales (typi- 
cally in the inertial range of scales) are described in a statistical consistent 
way (that is parameterized in terms of the resolved, large scale, velocity and 
temperature fields). This is done by filtering the governing equations for 
velocity and potential temperature by means of a filter operator. Applied, 
for example, to the ith-component of the velocity field, Ui, {ui = u, u 2 = v, 
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u 3 = w), the filter is defined by the convolution: 

h*w.-n* (4) 



UAX) 



where Ui is the filtered field and G(x) is a three-dimensional filter function. 
The field component Ui can be thus decomposed as 

Ui = Ui + u" (5) 

and similarly for the temperature field. In our model, the equation for the 
latter field is coupled to the Navier-Stokes equation via the Boussinesq term. 
Applying the filter operator both to the Navier-Stokes equation and to the 
equation for the potential temperature, and exploiting the decomposition 
(JHJ) (and the analogous for the temperature field) in the advection terms one 
obtains the corresponding filtered equations: 

dui dUiUj dr^f 1 1 dp 6 _ 2 _ 

-rrr = «— - ~TT a~ + 9i^~d i3 - fe ij3 Uj + vV m (6) 

Ot OXj OXj pOXi Oq 

g - ° _ < 7 > 

dt dxj dxj 

where p is the air density, p is the pressure, / is the Coriolis parameter, v is 
the molecular viscosity, k is the thermal molecular diffusivity, gi-^Sa is the 
buoyancy term and #o is a reference temperature profile. The quantities to 
be parametrized in terms of large scale fields are 



4 u) = UiU'j + u'luj + u'[u]- -rf = Qu] + Q"u 3 + Q"u], (9) 

that represent the subgrid scale (SGS) fluxes of momentum and heat, respec- 
tively. 
In our model: 

= -2K M {d i u j + d j u i )) (10) 

= -K H d t 6 (11) 

K M and K H being the SGS eddy coefficients for momentum and heat, re- 
spectively. 
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Table 1: The relevant parameters characterizing the simulated PBL. In this 
table, L x , L y and L z are the domain extension along the directions x, y and 
z, respectively; is the heat flux from the bottom boundary; U g is the 
geostrophic wind; Zi is the mixed layer depth, u* is the friction velocity and 
= Zi/u* is the turnover time; 



parameter value 



Lx-, Ly 


[km] 


2 




[km] 


1 


Q* 


[m K s" 1 ] 







[m s" 1 ] 


15 


Zi 


[m] 


461 


w* 


[ms" 1 ] 


0.7 




[s] 


674 



The above two eddy coefficients are related to the velocity scale e' , e' being 
the SGS turbulence energy the equation of which is solved in our LES model 
(Moeng, 1984), and to the length scale I = (AxAyAz) 1 ^ 3 (valid for neutrally 
stratified cases) Ax, Ay, and Az being the grid mesh spacing in x, y and z. 
Namely: 

K M = 0.117 1/2 (12) 

K H = 3^m- (13) 

Details on the LES model we used in our study can be found in Moeng, 1984 
and in Sullivan et al., 1994. Such a model has been widely used and tested 
to investigate basic research problems in the framework of boundary layer 
flows (see, for example, Antonelli et al., 2003 and Moeng and Sullivan, 1994 
among the others). 

3 The simulated PBL 

In order to obtain a stationary PBL we advanced in time our LES code 
for around six large-eddy turnover times, r*, with a spatial resolution of 128 3 
grid points. This time will be the starting point for the successive Lagrangian 
analysis (see next section). 

The relevant parameters characterizing our simulated PBL are listed in Table 
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Figure 1: The horizontally averaged velocity profiles, (a): stream- wise u- 
component, (b) span- wise t>-component. 



1 at t — 6 r*. At the same instant, we show in Fig. 1 the horizontally averaged 
vertical profile of the velocity components u, v. The average of the vertical 
component is not shown, the latter being very close to zero. We can observe 
the presence of a rather well mixed region which extends from z ~ 0.2 z^ to 
z ~ Zi. The energy spectra for the three velocity components are reported 
in Fig. 2. Dashed lines are relative to the Kolmogorov (K41) prediction 
E{k) oc A; -5 / 3 . Although the inertial range of scale appears quite narrow, 
data are compatible with the K41 prediction. 

4 Lagrangian simulations 

In order to investigate the statistics of pair dispersion, from the time t = 6r„ 
(corresponding to the PBL stationary state) we integrated, in parallel to the 
LES, the equation for the passive tracer trajectories defined by the equation 



We performed a single long run where the evolution of 20000 pairs has been 
followed starting from two different initial separations: -R(O) = Ax and 



x(t) = u(x(t),t). 



(14) 
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Figure 2: Energy spectra for the three components of the velocity field, (a): 
stream-wise, (b) span-wise, (c) vertical. The dashed lines correspond to the 
K41 prediction. 



R(0) = 2Ax, Ax being the grid mesh spacing whose value is 15.6 m. Tra- 
jectories have been integrated for a time of the order of 5000 s with a time 
step of around 1 s, the same used to advance in time the LES. 
At the initial time, pairs are uniformly distributed on a horizontal plane 
placed at the elevation Zi/2. Reflection has been assumed both at the cap- 
ping inversion (at the elevation Zi) and at the bottom boundary. 
For testing purposes, a second run (again started from t = 6 r*) with a 
smaller number of pairs (5000) has been performed. No significant differ- 
ences in the Lagrangian statistics have been however observed. The same 
conclusion has been obtained for a second test where the LES spatial resolu- 
tion has been lowered to 96 3 grid points. For a comparison see Figs. 3 and 4. 
The velocity field necessary to integrate (j!4j) has been obtained by a bilinear 
interpolation from the eight nearest grid points on which the velocity field 
produced by the LES is defined. 

In this preliminary investigation, we did not use any sub-grid model describ- 
ing the Lagrangian contribution arising from the motion on scales smaller 
than the grid mesh spacing. 
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t/x, 

Figure 3: The behavior of the (dimensionless) mean square relative disper- 
sion vs the (dimensionless) time. Full line: the initial separation is Ax; 
Dashed-line: the initial separation is 2Ax. Dotted line is relative to the t 3 
Richardson's law. 

4.1 Pair dispersion statistics 

In Fig. 3 we show the second moment of relative dispersion R 2 (t) for the two 
initial separations. Heavy dashed line represents the expected Richardson's 
law, which is however not compatible with our data for the largest initial 
separation 2 Ax. We can also notice how the R 2 (t) curve becomes flatter for 
larger separations. The same dependence has been observed by Boffetta and 
Celani (2000) for pair dispersion in two-dimensional turbulence. 

The fact that our data do not fit the Richardson law, for generic initial 
pair separations, is simply explained as a consequence of finite size effects 
(in space and in time) of our system. Indeed, it is clear that, unless t is 
large enough that all particle pairs have "forgotten" their initial conditions, 
the average will be biased. This is why we observe a consistent flattening of 
R 2 {t) at small times. Such regime is a crossover from initial conditions to 
the Richardson regime. From Fig. 3 we can see that the extension of such 
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crossover increases as the initial separation increases. 

Unfortunately, we cannot augment the time t too much because of the re- 
duced extension of our inertial range (see Fig. 2). To overcome this problem, 
and thus to allow a systematic estimation of the Richardson constant which 
does not depend on the choice of the initial pair separation, we use an alter- 
native approach based on statistics at fixed scale (Boffetta et a!., 2000). This 
is the subject of the next subsection. 

4.2 Fixed-scale statistics 

The characterization of transport properties in multi-scale systems, such as 
models of turbulent fluids, is a delicate task, especially when exponents of 
scaling laws and/or universal constants are to be measured from Lagrangian 
statistics. Additional difficulties arise in all cases where the standard asymp- 
totic quantities, for example the diffusion coefficients, cannot be computed 
correctly, for limitations due essentially to the finite size of the domain and to 
finite spatio-temporal resolution of the data. As we have seen in the previous 
subsection for the LES trajectories, the mean square relative dispersion, seen 
as a function of time, is generally affected by overlap effects between differ- 
ent regimes. We therefore use a mathematical tool known as Finite-Scale 
Lyapunov Exponent, briefly FSLE, a technique based on exit-time statistics 
at fixed scale of trajectory separation, formerly introduced in the framework 
of chaotic dynamical systems theory (for a review see Boffetta et al., 2000, 
and references therein). 

A dynamical system consists, basically, of a iV-dimensional state vector x, 
having a set of N observables as components evolving in the so-called phase 
space, and of a iV-dimensional evolution operator F, related by a first-order 
ordinary differential equations system: 



If F is nonlinear, the system (|15|) can have chaotic solutions, that is limited 
predictability, for which case an infinitesimally small error 5x on a trajectory 
x is exponentially amplified in time: 



with a (mean) growth rate A known as Maximum Lyapunov Exponent (MLE). 
The FSLE is based on the idea of characterizing the growth rate of a trajec- 



x(t) = F[x). 



(15) 



6x(t) ~ 5x(0) exp At 



(16) 
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tory perturbation in the whole range of scales from infinitesimal to macro- 
scopic sizes. In the Lagrangian description of fluid motion, the vector x is 
the tracer trajectory, the operator F is the velocity field, and the error 8x 
is the distance between two trajectories. It is therefore straightforward to 
consider the relative dispersion of Lagrangian trajectories as a problem of 
finite-error predictability. 

At this regard, the FSLE analysis has been applied in a number of re- 
cent works as diagnostics of transport properties in geophysical systems (see, 
for example, Lacorata et al., 2001; Joseph and Legras, 2002; LaCasce and 
Ohlmann, 2003). 

The procedure to define the FSLE is the following. Let r = \5x\ be 
the distance between two trajectories. Given a series of N spatial scales, or 
thresholds, 8%, 82, • • ■ ,8n have been properly chosen such that 8i+\ = p ■ 8i, 
for = 1, ■ • • , iV — 1 and with p > 1, the FSLE is defined as 

A <*> = WT) < 17 > 

where (T(8)) is the mean exit-time of r from the threshold 5 = 8i, in other 
words the mean time taken for r to grow from 5 to p5. The FSLE depends 
very weakly on p if p is chosen not much larger than 1. The factor p cannot 
be arbitrarily close to 1 because of finite-resolution problems and, on the 
other hand, must be kept sufficiently small in order to avoid contamination 
effects between different scales of motion. In our simulations we have fixed 
p = \/2. For infinitesimal 8, the FSLE coincides with the MLE. In general, 
for finite 8, the FSLE is expected to follow a power law of the type: 

X(8) ~ 8- 2h (18) 

where the value of 7 defines the dispersion regime at scale 8, for example: 
7 = 3 refers to Richardson diffusion within the turbulence inertial range; 
7 = 1 corresponds to standard diffusion, that is large-scale uncorrelated 
spreading of particles. These scaling laws can be explained by dimensional 
argument: if the scaling law of the relative dispersion in time is of the form 
r 2 (t) ~ t 7 , the inverse of time as function of space gives the corresponding 
scaling (JTSjl of the FSLE. In our case, indeed, we seek for a power law related 
to Richardson diffusion, inside the inertial range of the LES: 

\(8) = a8~ 2/3 (19) 
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where a is a constant depending on the details of the numerical experi- 
ment. The corresponding mean square relative separation is expected to 
follow Eq. ©• A formula can be derived, which relates the FSLE to the 
Richardson's constant (Boffetta and Sokoloff, 2002): 

a 3 / p 2/3 _ x \ 3 



p 2 / 3 ln p 



where (3 is a numerical coefficient equal to 1.75, e is the energy dissipation 
measured from the LES and a comes from the best fit of Eq. (|19|) to the 
data. Information about the existence of the inertial range is also given by a 
quantity related to the FSLE, the mean relative Lagrangian velocity at fixed 
scale that we indicate with 

v{8) = [<M5) 2 >] 1/2 (21) 

where 

5v(5) 2 = (£« - x^f (22) 

is the square (Lagrangian) velocity difference between two trajectories, x^ 
and x^ 2 \ on scale S, that is for la;* 1 ) — x^\ = 5. The quantity v(5)/5 is 
dimensionally equivalent to A(5), and, in conditions of sufficient isotropy, it 
represents the spectrum of the relative dispersion rate in real space. A scaling 
law of the type 

^ ~ 5-^ (23) 
o 

is compatible with the FSLE inside the inertial range and therefore with 
the expected behavior of the turbulent velocity difference as function of the 
scale. In Fig. 4(a) we can see, indeed, that the FSLE measured from the 
LES data follows the behavior of Eq. (|19|). from the scale of the spatial 
resolution to about the size of the domain. From the fit we extract the 
coefficient a = 0.1 m 2 ^t^ 1 . The energy dissipation measured from the LES 
is e = 6 • 10~ 4 m 2 t~ 3 . The formula of Eq. (J2*U|) gives a measure of the 
Richardson's constant C*2 ~ 0.5, affected, at most, by an estimated error of 
±0.2. In Fig. 4(b) we see, also, that v(5)/8 has been found very close to the 
behavior predicted by Eq. (J23)) . Variations within the error bars are observed 
by varying the spatial resolution from 128 3 grid points (triangles in Fig. 4) 
to 96 3 grid points (circles). 
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5/Lx 



5/Lx 



Figure 4: a) FSLE at two different resolutions. Triangles: 128 3 grid points; 
Circles: 96 3 grid points. The dashed line corresponds to a5~ 2 ^ 3 with a — 0.1 
m 2 t~ 3 . b) the same as in a) but for the relative velocity. The dashed line 
has slope —2/3. 
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5 Conclusions and perspectives 



We have investigated the problem of relative dispersion in a neutrally strat- 
ified planetary boundary layer simulated by means of Large-Eddy Simula- 
tions. In particular, our attention has been focused on the possible emergence 
of the celebrated Richardson's law ruling the separation in time of particle 
pairs. 

The difficulties in observing such behavior in a realistic PBL mainly rely on 
the fact that it is hard to obtain a PBL with a sufficiently extended inertial 
range of scales. For this reason, standard techniques to isolate the Richard- 
son's law and the relative constant turn out to be inconclusive, the results 
being strongly dependent, for instance, on the choice of the initial pair sepa- 
rations. To overcome this problem, we have applied, for the first time in the 
context of boundary layer physics, a recently established technique coming 
from the study of dynamical systems. As a result, a clean region of scaling 
showing the occurrence of the Richardson law has been observed and an ac- 
curate, systematic, measure of the Richardson constant became possible. Its 
value is C2 = (0.5 ± 0.2), where the error bar has been determined in a very 
conservative way. Such estimation is compatible with the one obtained from 
Fig. 3 in the case of initial pair separation equal to Ax. The important point 
is that the new strategy gives a result that, by construction, does not depend 
on the initial pair separations. As already emphasized this is not the case 
for the standard approach. 

Clearly, our study is not the end of the story. The following points appear 
to be worth investigating in a next future. 

The first point is related to the fact that in our simulations we did not use any 
sub-grid model for the unresolved Lagrangian motions. The main expected 
advantage of SGS Lagrangian parameterizations is to allow the choice of ini- 
tial pair separations smaller than the grid mesh spacing, a fact that would 
cause a reduction of the crossover from initial conditions to the genuine t 3 
law. The investigation of this important point is left for future research. 
Another point is related to the investigation of the probability density func- 
tion (pdf) of pair separation. In the present study, we have focused on the 
sole second moment of this pdf. There are, indeed, several solutions for the 
diffusion equation (JTJ) all giving pdfs compatible with the t 3 law. The so- 
lution for the pdf essentially depends on the choice for the eddy-diffusivity 
field, D(r). The answer to this question concerns applicative studies related, 
for example, to pollutant dispersion because of the importance of correctly 
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describing the occurrence of extreme, potentially dangerous, events. 
Finally, it is also interesting to investigate whether or not the Richardson law 
rules the behavior of pair separations also in buoyancy-dominated boundary 
layers. In this case, the role of buoyancy could modify the expression for 
the eddy-diffusivity field, D(r), thus giving rise to an essentially new regime 
which is however up to now totally unexplored. 
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